close all;
clear all;

%% Please modify this part
data_root = 'D:\2011\LPCVD_2010.08.27\Simulations\2DSPDE_EW\closed_loop\2DEW_ModelFitting4_CL\XYmax_8000nm_Stronger_Constraints\r2_and_m2\';
data_subfolder = {'qr_0.00000001-qm_1.0';
                  'qr_0.0000001-qm_1.0';
                  'qr_0.000001-qm_1.0';
                  'qr_0.00001-qm_1.0';
                  'qr_0.0001-qm_1.0';
                  'qr_0.001-qm_1.0';
                  'qr_0.01-qm_1.0';
                  'qr_0.1-qm_1.0';
                  'qr_1-qm_1.0'};
weighting_ratio = [0.00000001,0.0000001,0.000001,0.00001,0.0001,0.001,0.01,0.1,1];

%% 
r2_end = zeros(size(weighting_ratio));
m2_end = zeros(size(weighting_ratio));

for i = 1:length(weighting_ratio)
    record_r2 = SinkToFile([data_root,data_subfolder{i},'/Output_r2_nm2_1.dat'],1);
    record_r2.VariableName = '<r^2>';
    record_r2.VariableUnit = 'layer^2';
    
    record_m2 = SinkToFile([data_root,data_subfolder{i},'/Output_m2_1.dat'],1);
    record_m2.VariableName = '<m^2>';
    record_m2.VariableUnit = '';
    
    r2 = record_r2.mean;
    r2_end(i) = r2(end);
    m2 = record_m2.mean;
    m2_end(i) = m2(end);   
end


[AX,H1,H2] = plotyy(weighting_ratio,m2_end,weighting_ratio,r2_end);
set(AX(1),'position',[0.18 0.15 0.675 0.775]);
set(AX(2),'position',[0.18 0.15 0.675 0.775]);
set(AX(1),'Xscale','log');
set(AX(2),'Xscale','log');
set(AX(1),'FontSize',20);
set(AX(2),'FontSize',20);
set(get(AX(1),'Ylabel'),'String','<m^2_{\Delta}>_{ss}','FontSize',20);
set(get(AX(2),'Ylabel'),'String','<r^2_{\Delta}>_{ss}(nm^2)','FontSize',20);
set(get(AX(1),'Xlabel'),'String','q_{r^2}/q_{m^2}','FontSize',20);
set(H1,'LineWidth',2,'Marker','o');
set(H2,'LineWidth',2,'Marker','s');
legend(AX(1),'m^2','location','nw');
legend(AX(2),'r^2','location','se');
print('-depsc2','r2_m2_vs_weighting.eps');

figure();
subplot(2,1,1);
plot(weighting_ratio,m2_end,'LineWidth',2,'Marker','o');
set(gca,'Xscale','log');
set(gca,'FontSize',20);
ylabel('<m^2_{\Delta}(t_f)>','FontSize',20);
xlabel('q_{r^2}/q_{m^2}','FontSize',20);
%legend('m^2','location','ne');

subplot(2,1,2);
plot(weighting_ratio,r2_end,'LineWidth',2,'Marker','s');
set(gca,'Xscale','log');
set(gca,'FontSize',20);
ylabel('<r^2_{\Delta}(t_f)>(nm^2)','FontSize',20);
xlabel('q_{r^2}/q_{m^2}','FontSize',20);
%legend('r^2','location','ne');
print('-depsc2','r2_m2_vs_weighting_v2.eps');

%% Contour plot
r = sqrt(r2_end)*1e-9;
m = sqrt(m2_end);
m_mesh = linspace(1*min(m),1.1*max(m),30);
r_mesh = linspace(1*min(r),1.1*max(r),30);
R_mesh = zeros(length(r_mesh),length(m_mesh));
for i=1:length(r_mesh)
    for j = 1:length(m_mesh)
        R_mesh(i,j) = Reflectance(r_mesh(i),m_mesh(j));
    end
end

%%
figure();
contourf(m_mesh,r_mesh*1e9,R_mesh);hold on;
plot(m,r*1e9,'ro','MarkerSize',8,'MarkerFace','w','MarkerEdge','k');
colorbar();
box on;
set(gca,'FontSize',20);
xlabel('m_{\Delta}');
ylabel('r_{\Delta} (layer)');
print('-djpeg','r2_m2_simultaneous_contour');   
print('-depsc2','r2_m2_simultaneous_contour.eps');   






